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We provide a comprehensive theoretical framework and a quantitative test of the 
method we recently proposed for processing data from a spherical detector with 
five or six transducers. Our algorithm is a trigger event generator performing a 
coherent analysis of the sphere channels. In order to test our pipeline we first built 
a detailed numerical model of the detector, including deviations from the ideal case 
such as quadrupole modes splitting, and non-identical transducer readout chains. 
This model, coupled with a Gaussian noise generator, has then been used to produce 
six time series, corresponding to the outputs of the six transducers attached to the 
sphere. We finally injected gravitational wave burst signals into the data stream, 
as well as bursts of non-gravitational origin in order to mimic the presence of non- 
Gaussian noise, and then processed the mock data. We report quantitative results 
for the detection efficiency versus false alarm rate and for the affordability of the 
reconstruction of the direction of arrival. In particular, the combination of the two 
direction reconstruction methods can reduce by a factor of 10 the number false alarms 
due to the non-Gaussian noise. 

PACS numbers: 04.80.Nn,95.55.Ym 



I. INTRODUCTION 

Started more than forty years ago with the pioneering work of Weber [Tj, the rush for 
direct gravitational wave (GW) detection have nowadays reached a crucial phase. On one 
side, the era of large interferometers (LIGO, VIRGO, GEO, TAMA) [2] has begun and the 
next, so-called "advanced" versions of the interferometers are currently considered as the 
most likely candidate to make the first detection (although uncertainties in stellar population 
estimates still do not allow to keep this goal for granted, see 
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On the other side, resonant bars [3] hke the INFN ones (EXPLORER, NAUTILUS, 
AURIGA) or ALLEGRO seem to have exhausted their original leading role after several 
years of impressive performances (see for instance [3J for some reports on the most recent 
runs) which made them by far the best instruments for a long time. 

This does not mean at all that resonant detectors are going to be out of business, as 
resonant bars are currently running as a network [6j and joint analyses with interferometers 
have already taken place [7j (and are foreseen in the future). Moreover high-performance, 
new generation antennas (like large spheres [S] or dual detectors [9]) have already been 
conceived. In this framework, the two small spherical detectors in Leiden [TU] and San 
Paulo pjj, that are already in the commissioning phase, can rightfully be considered as the 
first specimens of this new generation. 

What makes a sphere really different from bars and interferometers is the fact of being 
a multichannel detector. This distinguishing feature not only determines its isotropic sen- 
sitivity and its capability to reconstruct the GW direction, but also opens the way to new 
opportunities, and issues, from the data analysis point of view. After the seminal work of 
Wagoner and Paik [12], where the basic features of the detector have been pointed out, more 
and more detailed sphere models and configurations have been studied by several authors 
[13] [TU [T5l [T6l [T8l \T9\ [20| [2T| [22] . In particular, various solutions have been proposed to 
the problems of the ideal transducers configuration |20], parameter reconstruction |13l [T9] . 
and multidimensional data analysis [TU] . 

The analysis method that we have set up and studied in the present work is the result 
of a synthesis and a refinement of some of these contributions, and is intended to represent 
the core of the pipeline that we are building for the miniGRAIL detector. Some preliminary 
results have already been proposed in [IT], here we expose the theoretical foundations and 
assess quantitatively the solidity of the method as an event trigger-generator. 

Particular care has been dedicated to all those concrete details that enter the game as 
"theoretical" data analysis is turned into a real, working device: these include, in primis, 
attention to the problem of false detections and keeping under control the computation time 
required by the pipeline. 

The first part of the paper is devoted to building an accurate detector model and to 
the generation of a set of simulated data, expanding the treatment in PT]. We consider 
a real, imperfect sphere with six read-out electromechanical amplificators of its surface 
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oscillations. Our method works also with five outputs, while using a smaller number of 
transducers significantly degrade the performances of the detector, even if four transducers 
are in principle sufficient to fully reconstruct the GW tensor. 

After having written the basic equations, we work out the spectrum of each transducer's 
output and we generate simulated data starting from the elementary noises which enter the 
detector at different levels. Then we discuss the problem of extracting the GW parameters 
out of the data. 

In the second part we describe the method to generate triggers of events, based on the 
coherent WaveBurst algorithm, the event trigger generator used by the LIGO- VIRGO burst 
data analysis group. We adapted it to a multimode detector, where different channels 
have correlated noise, differently from the LIGO- VIRGO situation where the detectors are 
not co-located. The availability of several channels allows a determinaton of the arrival 
direction of the burst through a standard likelihood analysis method. Such a determination 
is then validated by a consistency check with the geometrical method exposed in paragraph 



III C , named determinant method, which provides an independent indication of the arrival 



direction. 

In the last section the efficiency of the method for four different amplitudes of injected 
bursts is shown. We tested the method against injections of both GW and non-GW signals, 
where for non-GW signals the injection amplitude is shared randomly the five quadrupolar 
channels. 

While each of the two above-mentioned method alone is not able to distinguish between 
GW and non-GW excitations, we find that the combination ot the two allows to reduce the 
false alarm rate by roughly one order of magnitude, while keeping a good detection efficiency, 
and provides a good determination of the arrival direction of the signal. The multi-mode 
coherent analysis is thus able to compute all the relevant parameters of a GW and to reduce 
the false alarm rate. 



Figure 1: The readout scheme, from (TTj. 



II. MODEL OF THE DETECTOR 
A. Equations for modes, transducers and readout currents 

The sphere is weU modeled as a set of coupled oscillators, which describe the dynamics of 
the relevant sphere vibrational modes, of the transducers, and of the electrical circuit that 
are at the core of the readout devices. As the relevant equations have been already discussed 
in detail by several authors (see for instance [IHl 1221 ES]), we jump to the mathematical core 
of the problem, skipping introductory material and definitions that can be found in the 
literature. These equations hold for any number of sphere modes and transducers, although 
in the following subsections we will specify our choices for such numbers. For the sake of 
clarity and completeness, the definitions and values of all the quantities appearing below are 
summarized in figjl] and in Appendix A. 

In Fourier space, the amplitudes C,n of the various radial modes of the sphere obey to the 
following equations 



2 .WfcCjX rrik . Ek 



In - aN ^Nk^Ji 



(1) 



,^NkJ^I 

k 

which relate their dynamics to that of the transducer displacements qk and of the primary 
currents through the position matrix i?jvfc = ^Ar(6'fc, 0^)^32]. The equations for the trans- 
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ducer displacements are: 



9 V — V / 9 9 \ -^k i / \ 

-U y BkNOtNiN + \^^k- ^ +i-p^](lk-i Jk = fk, (2) 

^ V Qk J rukuj 

while the readout circuit can be described through the equations driving the primary currents 
and the secondary ones, I^, which are proportional to the actual outputs of the detector: 



EkQk + 



J,-zuMl\h-fl) = fl, (3) 



- icuMf J, + {LI + LI) (4 - fl) = n ■ (4) 

The /'s appearing above are the stochastic forces related to the various dissipative com- 
ponents of the detector [33]; their action can be described through the corresponding noise 
spectral densities: 

< Uu)f;{u;') >^ S,yiu)5iu - iu') . (5) 
The non- vanishing components of the noise spectral density matrix are: 

Snn' = '^^bT—j^—Snn' , Slf., = 2kBT ^ 5kk' , S^^, = 2kBTrkSkk' , (6) 

MC^jV ^k^k 

and, according to the Clarke model [21] for the SQUID, 



Slk'i^) = llu;'^{M',^)%,, , (7) 
^k 

where Snn', Sl/^,, S^f^,, S^j^,, Sl/^,, S^'l, are the spectral densities of, respectivley: the forces 
acting on the the radial mode displacements, on the transducers, on the primary currents, 
on the secondary currents, the spectral densities of the secondary current noises and the 
spectral density of the secondary current noise forces times the secondary current noise. 
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B. Output current determination and simulation 

1. Specification of the system 



By using the following matrix notation 



Z = 





-V [ajv] • BNk ■ 


V [ajv] • BNk ■ [sk/ruk] 
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(8) 



(9) 



Q = {^N: Qk: Jk, h - fk) : ^ = {/iV, /fe, /fe , fk} : 



(10) 



where X„ is the n-dimensional identity matrix, the dynamical system of equations can be 
written in a concise way 



(11) 



which in turn implies 



{Cat, Qk, Jk, Ik}^2-'-A-J^+ {0, 0, 0, fl} . 



(12) 



This reduces the problem of finding the currents Ik, and their spectral correlation matri- 
ces, to the inversion of the matrix Z. 

In some particular case, such inversion is very easy to compute. For instance, for the 
case of 6 transducers in the TIG A configuration, if we include in the model only the first 5 
quadrupolar sphere modes and the scalar one, then 
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BNk ■ - BIj^ ■ Bnu - ^1<3 ■ (13) 

In this case, in the ideahzed hypothesis that all the 6 transducers and readout chains are 
identical, thus implying that all the blocks of the kind T) [X^] appearing in Z are proportional 
to the identity, one obtains that the following redefinition 



with 



leads to 



7^ 
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= n 


■Q', 




It 

-Ln 














BNk 














BNk 
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BNk 


Z-Q 


= 7^ 


■Z'- 


Q', 



(14) 



(15) 



(16) 



where all the 16, 6 x 6 blocks of Z' are diagonal. 



Since TZ is easily invertible by means of eq.(13), the J^'s are easily found. This particular 



case is nothing but a rephrasing of the idea of mode channels of an ideal TIGA configuration, 
see [IH][31]. 

However, since our purpose is to simulate a realistic detector, we will not take advantage 
of any idealization, and we will invert the Z matrix numerically. 

For the same reason we have decided to include many sphere modes in our model, con- 
trarily to the common practice of considering only the lowest quadrupole multiplet, that is 
the one that interacts with the GW. More precisely, we have included in our model all the 
radial modes with resonant frequency below that of the scalar mode: this gives a total of 
30 degrees of freedom in this part of the system, distributed among one scalar, one vector, 
two quadrupolar, one 1=3 and one 1=4 multiplets. The precise reason of this choice will be 



discussed in subsection II C for the moment we just mention that even if the extra modes 
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Figure 2: Theoretical output of one of the transducers when 6 or 30 sphere vibrational modes are 
included in the model. To emphasize the differencies, a restricted region of the spectrum have been 
displayed: this is the reason why only nine out of the eleven resonant modes are visible. 

resonate far away the interesting region of the spectrum (which is around 3kHz), neverthe- 
less their inclusion in the model affects the response function even in the relevant spectral 
region, shifting the peaks in the response function of a few Hz's, as can be seen in figj2j In 
a real detector, however, such shift may be masked by the action of the suspension system, 
finite-size effects of the transducers and/or some other unmodelled factors. At the moment 
we are not able to predict the magnitude of these unmodelled effects for miniGRAIL. 

For what concerns the transducers arrangement, we will stick to the TIG A configuration 
(with non identical transducers, however), although our analysis can easily be adapted to 
other configurations, provided the number of transducers is kept equal to six. 

2. From expectation values to a typical noise realization 



Once Z has been numerically inverted, we can proceed to determine the currents Ik through 
eq.([l2|, and in particular we can estimate the corresponding spectral density matrix from 
eqns.(§-(|7|): 

< h{oo)I*Au;') SUu;)6{u; - u') . (17) 
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Figure 3: The calculated output of the transducer ^^0, along with the various noise contributions. 
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Figure 4: Simulated output of transducer ^^0, from [T7] . 



The square root of one of the diagonal elements of 5*^^, is shown in figjsj along with the 
different noise contributions. 

This is still not enough for our purposes: what we really want to get at this stage is 
not just expectation values and cross-correlations, but an actual realization of a possible 
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detector output. To obtain this, we have generated several time series of white noise, we 
have taken their numerical Fourier transforms, and " colored" them according to the spectral 
dependencies of eqns.([6])-([7]), thus obtaining a possible realization of the Fourier transform 
of the various stochastic forces /'s. Finally, we have derived, still in Fourier space, the 



realizations of the currents Jfc, by means of eq.(12), satisfying eq.(17); one result of this 
procedure is displayed in figjlj 

We underline the fact that fig|4] does not represent just a white spectrum "colored" 
according to the expectation value shown in fig|3j but it rather contains, in virtue of the 
procedure we have adopted to generate it, all the correct cross-correlations with the other 



output currents, i.e. the non-diagonal elements of eq.(17). This fact is crucial if the injected 
signals are to be reconstructed. 

To summarize the results of this subsection, we have numerically generated the six current 
outputs of a realistic sphere working in the TIGA configuration. Our model is realistic in 
the sense that we have chosen the values of all the components of the simulated detector 
according to what is actually being implemented on miniGRAIL, and we have included an 
adequate number of sphere vibrational modes. 



C. Interaction with the GW: the transfer matrix 



We now change point of view: we assume that the six secondary currents are given and 
we would like to determine the sensitivity of our detector to GW's. In other words, we are 
trying to reproduce as faithfully as possible a real experimental situation. 

We remind that when a GW impinges the detector, a deterministic force is added to the 
quadrupolar components of /at: 

fN^ fN + ^Rx'hN for Ar = 0,...,4, (18) 

where R is the sphere radius, x = —0.3278 for CuA16% (the alloy miniGRAIL is made 
of) and the h^^s are the quadrupolar spheroidal components of the GW, defined as follows 
through decomposition via the tensor spherical harmonics: 



hi 



Af=0,4 



3^^'^ --h 



N 



(19) 
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with niN = 0, Ic, Is, 2c, 2s and Yml = 'Yli jy^i^,ij^^^^ i being ra* the versor of the arrival 
direction of the wave. 



1. Issues in getting the h^'s from the transducer outputs 

The problem is clearly to get the most possible information about the /iat's from the knowl- 
edge of the six /fc's. For this goal, it would be enough to determine, and most of all to invert 
somehow, the T^n matrix, which is defined as the k x 5 block determined by crossing the 



last k lines with the first 5 columns of the matrix Z which appears in eq.(12). Basically, 
TkN is the matrix entering the relation: 

h = TkN ■ In ■ (20) 



with fisf given by eq.(18), which holds in the limit of very large SNR. 

In general, the solution to this problem is very much dependent on the number of trans- 
ducers. If such number is less than four, then there are not enough experimental inputs to 
determine all the /iat's, and thus the matrix hij cannot be completely reconstructed. Rather, 
only some linear combinations of the h^'s (or, in another notation, of the GW polarization 
components) can be determined |35j. 

If there are four transducers, then in principle the metric perturbation can be recon- 
structed, provided one imposes a priori on the data the constraint that the five h^^s should 
form a transverse tensor. The five transducers case is similar, with the advantage that one 
can use the transversality condition as a veto, thus reducing the false alarm rate at a given 
SNR, rather then imposing it from the beginning. 

If the transducers are six, as in the present case, the situation is more involved because 
the system is over-constrained by the fact that the 5 /ia/'s are to be determined by the 6 
outputs. This kind of problem has generally no solution, as it is physically correct since the 
Jfc's are not determined only by the GW, but also by the detector noise; however in some 
cases an optimal solution can be found. 

First of all, we have seen that in an ideal TIGA configuration (identical readout chains) 
a change of basis brings the matrix Z' to a block diagonal form. By inspection of the form 



of the system of equations (11) it can be seen that this implies also that the T^n matrix has 
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kN 



(21) 



the following pseudo-diagonal form: 

/* o\ 
0*000 
0*00 
*0 
* 
\0 / 

the last row of zeroes meaning that the sixth redefined component of is independent of 
the /iat's. In this case one has just to consider the square, and invertible, matrix made of 
the first 5 lines of TkN'- i-e. only the first five J^'s (which are nothing else but the TIGA 
mode channels) are used to find the five /iat's and the system is again well determined. 

When the idealized conditions for a perfect TIGA configuration no longer hold, the 
situation becomes more complicated because there is no simple way to get to the form 



(21) for the T^tv matrix; in particular, in the general case all the six current outputs (or any 
linear combination of them) depend on the GW parameters. 



In this case, the system (20 ) is really over-constrained and has in general no exact solution; 



in such cases, the best one can do from the mathematical point of view is to find the best 
possible approximate solution for the system, that is, in the case at hand, a vector such 
that the norm of the vector 



N 5 



(22) 



is minimal. If the matrix T^^ has maximum rank, then the pseudo solution f"^ can be shown 
to be unique, and to be given by 



(23) 



This strategy, which has been proposed for example in [19] is certainly a good one in the 
case of very high SNR, while at low SNR, and in particular in the determination of the 
detector's sensitivity (where SNR = 1 by definition), it shows some drawbacks. 

The problem is that, when the noise gives a non negligible contribution to the output 
currents, there is no physical reason to expect that the actual /iat's are the ones that minimize 



the norm of the vector (22) 
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2. Inclusion of the scalar mode 



Our proposal in this case is to include the scalar vibrational mode of the sphere, so that 



eqns.(18), (19) are now turned into 

1 



fN ^ fN + ^RXNu:''hN for iV = 0,...,5, (24) 



h^^= E 3^S..^^ + #^5, (25) 



N=OA 



being Xjj the identity matrix, and xn = X for = 0, ... ,4, while it can be shown that 
X5 = -3.8. 

As we have already discussed, inclusion of the scalar mode makes the position matrix B^k 



invertible and orthogonal, see eq.(13), because the row corresponding to the scalar mode is 
orthogonal to the 5 lines related to the quadrupolar modes. This means that adding the 
scalar mode is the best possible way to complete the 5 quadrupole position vectors to a basis 
of the six-dimensional linear space. 

If we proceed this way, another sixth column must be included in the definition of 7/vfc, 
which now becomes an invertible 6x6 matrix. In other words, including the scalar mode is 
a way to turn an over-constrained rectangular system into a squared one; it is clear that the 
same result could have been obtained by adding any other sphere vibrational mode rather 
than the scalar, but as previously stated, inclusion of the scalar mode makes the system 
particularly symmetric. 

A word of caution is needed here: by the inclusion of the scalar mode, and in particular 



its parameterization by eq. (25), we do not aim at detecting scalar gravitational radiation. 
We are rather parametrising the noise contained in the scalar vibrational mode of the sphere, 
and dropping it from the analysis of the quadrupolar modes: in other words it makes much 
sense to take the part of noise related to the trace out of the 5 quadrupolar /iat's, and the 
addition of the scalar mode, which is sensitive exactly to this part, is the way to do it. 



Stated in yet another way: the system (20) is over-constrained, and one needs an extra 



input to obtain the best possible solution for h^. In the idealized TIGA case, this extra 
input comes from the existence of mode channels, while in the high SNR case a good criterion 
is the minimal distance one. In general, we think that a physically motivated option is to 
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express the six noisy output currents in terms of the five quadrupolar sphere modes plus 
the scalar one, that is the best way to complete the quadrupolar modes to a basis in the 
six-dimensional space, and to ensure that the trace part of the noise does not "leak" into 
the quadrupole modes. 



3. Transfer matrix 



We now come to the determination of the transfer function and of the GW components. 



The system ( 26 ) can be rewritten as follows 

h{u:) = ]^Ru:^%N{uj)-V[xN]-hN{u:), N,k = 0,...,5, (26) 

and numerically inverted (for every value of uj) to give 

hN^uj) = TNk{^^) ■ Ik{uj) , (27) 

where T/vfc(t^) defines the transfer function (which in this case is actually a transfer matrix) 
of the system. 

Eq.(27) can be combined with the results of the previous section, in order to give the 



spectral density matrix, which can be used to estimate the detector sensitivity. The left part 
of fig. [5] shows the diagonal components of such sensitivity matrix, that is square root of the 
expectation values < h^h*^ >, for all values of A^. Since the different modes have different 
sensitivity curves, and since they are also cross-correlated, the overall sensitivity of the sphere 
is not isotropic. On the right part of the same figure we have shown such sensitivity averaged 
over all the incoming directions and over all the possible linear polarizations, together with 
the various noise contributions, and we have compared it with the corresponding quantity 
for the LIGO interferometer [36] . 

It can be noticed that i) the five quadrupolar components having different sensitivities is 
due to the fact that the readout chains are not identical and that the mode resonances are 
not degenerate, and ii) as expected the detector is poorly sensitive to the scalar mode: this 
happens because the transducers are tuned to the quadrupolar multiplet whereas the scalar 
mode resonance is at a much higher frequency then the quadrupole (roughly at 5kHz). 
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Figure 5: On the left, all the GW quadrupolar modes and the scalar one. On the right, the sky- 
averaged sensitivity of the sphere (along with the various noise contributions), compared to that 
of LIGO. 



Coming from expectation values to our detector simulation, the same procedure can be 
applied to the Fourier transform of the Jfc's realization produced in the previous subsection, 
to give the Fourier transform of the h modes. 

For example, the quadrupolar /iq mode spectrum is shown in fig. [6} together with its 
expectation value. 

Once the hi^{uj)^s are transformed back in the time domain, six time series are obtained, 
which describe the noise corresponding to the components of the 5 spheroidal modes of the 



tensor hij and of the scalar mode, see eq. (25). However, given the poor sensitivity of the 



detector to the latter (see fig. |5]), this channel is bound to contain instrumental noise only 
and it will be dropped in the forthcoming analysis. 



III. THE ANALYSIS PIPELINE 

We have now the five time series corresponding to the quadrupolar modes, which are 
enough and necessary to reconstruct the most general symmetric, traceless 3x3 matrix. 

Still this is a redundant description of a GW, and we will exploit this redundancy to 
discriminate between real GW signals and excitations of the modes due to disturbances 
other than gravitational. 

The treatment in this section closely follows the exposition of [17], which is summarized 
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Figure 6: The mode /iq as it has been reconstructed starting from an actual noise reahzation, 
superimposed with its expectation value. 



here. 



A. The scalar trigger 

If one knew exactly the form of the expected signal, the optimal strategy would be to 
perform a multidimensional matched filter as follows: 



(28) 



where the /i^^(a;)'s are the Fourier modes of the quadrupolar components of the expected 
signal, and we have introduced the following notation: 



(X ■ Y)s = Xl,ico) ■ iS%'^,iu) -Y^icu) iV, AT' = , . . . , 4 . 



(29) 



However in a realistic case neither the temporal dependence of the signal, nor the arrival 
direction are known. One could think of overcoming this difficulty by doing several matched 
filters with signal of various shapes and polarizations, and coming from several directions 
in the sky, but this method is computationally very intensive, and in addition results in a 
much higher false alarm rate than the method we are going to describe now. 
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We would like to find triggers out of the stretch of data, without need of a detailed 
knowledge of the signal. This has been done in [T7j by building the following quantity: 

Hiuj) = {h-h)s. (30) 

The Fourier transform of this quantity is then fed to a suitably adapted version of Wave- 
Burst, the event trigger generator used in the LIGO data analysis The WaveBurst 
algorithm make use of the wavelet decomposition, and among the bank of wavelet packets 
we picked the Symlet base with filter length sixty [31j. 

If a GW signal is present in the h^'s, WaveBurst will detect it in H whatever its shape, 
polarization and arrival direction, provided of course, that it is strong enough. 

Once the trigger has been established the analysis is performed on the modes /^at's, by 
collecting for each trigger the values of the wavelet coefficients. At this point the arrival 
direction can be reconstructed by a likelihood method analysis. This is different than what 
is being done by the Coherent WaveBurst algorithm [30], where the likelihood value is 
computed at every point of the time-frequency plane, as we do not have the sufficient 
computational power to perform such a daunting task. On the other hand, as it will be 



shown in sec. IV, adopting the above mentioned scalar trigger allows a good detection 
efficiency and reduce enormously the computational cost of the procedure, as our algorithm 
runs on a standard MacPro machine with two 3GHz processors and it takes for the analysis 
a time which is roughly two thirds of the actual time duration of the data. 

The values of ^, of the GW arrival direction are then found as the ones maximizing the 
likelihood function. This method identifies an arrival direction of the candidate GW event 
no matter if a real GW has excited the detector or a glitch, say, has taken place. To confirm 
this direction determination we combine it with a different, algebraic method. When the 
two methods do not determine the same directions, within some tolerance to be discussed 



quantitatively in sec. IV, the event can be discarded as spurious. 



B. Likelihood method for direction reconstruction 



The probability of having a given stretch of data {/lAf} assuming a GW with polarization 
amplitudes /i+ and hx from the direction identified by the usual polar angles 6, 0, is often 
referred as likelihood function and is denoted by p{{h]\f}\H(^j^). For stationary, Gaussian, 
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Figure 7: Example of event trigger. Sine-gaussian GW injection, centered at t = 500sec, width 
50msec, central frequency 2825Hz and strength hrss = 7 x lO^^^Hz^^/^. On the left is the H mode 



defined in eq. (30), and on the right the likelihood of eq. (34) on the points of the trigger. 



white noise with zero mean, and by taking into account that the noises in the different 
channels are correlated, such probability is 

p{{hr,}\H^,) (X exp [-{{h -0-{h- 0)5/2] , (31) 

where Ctv is defined as 

CN = F+{9,^)h+ + F^{9,(P)h^, (32) 

Fp^ being the pattern function of the mode for the +, x polarization. 

Following the standard procedure [271 128] , the likelihood ratio A can be defined as 



p{{hN}\Ho) ■ 



(33) 



By maximizing the likelihood ratio eq. (33) with respect to h^, a function C{9, (j)) of the 
angles alone is found 



1 



[(F+-F+)5(FX-FX)5-(F+-FX)|] 



21-1 



X 



(34) 



[FX ■ F-)s{h . F+)| + (F+ ■ F+)s{h - FX)! - 2{h ■ F+)s{h ■ F-)s{F+ ■ F^)s\ . 

Maximizing the likelihood ratio is equivalent to maximizing the posteriori probability 
p{H(^j^\{hN}) in the case of fiat priors. We then take the values of 6', (which enter the 
expression for the F/v's) maximizing Xlitriooerl 



{trigger} ^h\0 , (p), 



that is the sum of the likelihood 



over every point of time-frequency plane exceeding the threshold, as our estimation of the 
arrival direction of the candidate event. I 
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C. Determinant method for direction reconstruction 

Another method to reconstruct the direction of arrival GW by exploiting the (redun- 
dancy of the) five quadrupolar modes is based on linear algebra considerations [18J. A 
general metric perturbation in the transverce-traceless gauge is necessarily parametrized by 
a symmetric, traceless 3x3 matrix with zero determinant. No matter which arrival direction 
nor polarization, the null eigenvector always points to the arrival direction (with the usual 
ambiguity in up and down direction). 

As the quadrupolar modes allow full reconstruction of the metric perturbation, we have 
an independent method to identify the incoming direction of the signal. 
Of course none of the eigenvalues is expected to be exactly zero at every point in the trigger, 
we then proceed to order the three eigenvalues Aj so that |Ao| < |Ai| < IA2I and define the 
following quantity 



r = 



(35) 



v/AfTAf 

For a perfect GW-like r vanishes, thus the smallest it is, the less the noise is contaminating 
the GW signal. If for at least one point in the trigger r < tq, with tq a suitably chosen 
threshold, a direction can be identified by associating to the signal the direction designated 
by the eigenvector relative to the smallest eigenvalue Aq. 

For each point of the trigger satisfying the condition r < tq a direction is identified; the 
average is then taken by weighting them with the factor 1/r in order to obtain a single 
direction for each trigger. 



IV. TEST OF THE METHOD 
A. Setup and calibration 

In order to test the efficiency of our method, we injected mock GW signals in the five 



channels corresponding to the spheroidal components of hij (see subsection II C). 

As efficiency is a good measure of the validity of the method only in presence of an 
estimate of the false alarm rate, we also injected non-GW signals which are supposed to 
mimic the presence of non-Gaussian noise in the data. 
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Figure 8: The stars indicate the injections, while the dashed line is the strain sensitivity for 
the arrival direction of the injections. For comparison, we have displayed also also the averaged 
sensitivity (continuous line), as well as the ones corresponding, at any given value of the frequency, 
to the best and worst possible directions. 



In both cases the amplitude shape is a sine Gaussian 

h,^,{t + to) = hoe-'" /^^''"hm{2n ft) 



(36) 



where A = O.OSsec, / = 2825Hz, to is the signal center and ho is determined by the signal 
hfgg. 

Four sets of injections (of 50 signals each) have been done both for GW and non-GW sig- 
nals, with amplitudes hrss = {3, 5, 7, 10}-10^^^Hz~^/^, corresponding to SNR ^ 16, 27, 37, 53 
in amplitude, see figure [8j 

For the mock GW signals, the amplitude distribution among the different channels corre- 
sponds to linearly polarized GW's (either + or x) with arrival direction 6 = 54.736°, = 45°, 
while for non-GW signals the amplitude is distributed randomly on the different channels 
and normalized to ensure the required hrss- 

For each injection detected by the trigger, the arrival direction has been reconstructed 



by the two methods (likelihood and determinant) explained in sec. Ill and the event have 



been selected only when the two methods had a discrepancy below some threshold 6. 

A ROC curve can be constructed to assess what is the best value for 6: for each set 



21 



of injections one can define tlie efficiency as the fractional number of injections which are 
recovered and for which the two different direction identification methods give directions 
separated by less than 6. By measuring the efficiency for GW and non-GW injections and 
plotting the former vs the latter for different values of 6, fig. [9] has been obtained, where the 
four curves refers to the four different signal strengths. 

The same can be done comparing the GW detection efficiency with the false alarm rate 
obtained analyzing the data containing only Gaussian noise (i.e. without any injection): the 
plot is reported in fig. 10 and shows that, contrarily to what happens in fig. [9| the ROC 



curve shows a strong dependence on 6, especially at low hrss (notice that the scales of the 
two figures are different). 

This calibration enabled us to set at 6thr = 0.2 rad as the optimal value; all the results 
reported in the remaining part of the paper have been obtained with 6 = 6thr- 

B. Results 

1. Likelihood only 



As a first step, we want to check what is the significance of finding a trigger on the scalar 



mode H defined in eq. 30 In Fig. 11 we report the distribution of the likelihood values for 
the triggers obtained in correspondence of GW and non-GW injections. Mean values and 
sigmas for the likelihood are reported in tab. [Tj 

We see that our algorithm detects more efficiently GW injections rather than non-GW 
ones, but at this stage it is impossible to discriminate a priori between the two kinds of 
injections. 

2. Likelihood + determinant 



We now complete the maximum likelihood method with a cross check in order to get rid of 
spurious events: this is achieved by means of the geometrical method described in subsection 

[mcl 
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Figure 9: Detection efficiency of GW-injections vs. detection efficiency of non-GW signals ob- 
tained by varying the maximum allowed distance between the two different direction reconstruc- 
tion methods. Lines of increasing thickness correspond to injections sets of increasing amplitude, 
hrss = {3,5,7,10} X 10-2iHz-^/2 




False alarm rate (Hz) 



Figure 10: Detection efficiency of GW-injections vs. false alarm rate of events on the data with 
just intrinsic detector noise, obtained by varying the maximum allowed distance between the two 
different direction reconstruction method. Lines of increasing thickness correspond to injections 
sets of increasing amphtude, hrss = {3,5, 7, 10} x 10~^-'-Hz~^/^. 
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Figure 11: Distribution of the values of the likelihood function for GW-injections (empty-black) 
and for non-GW injections (filled-green) for various h^ss- 
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a 


10-20 


50/50 


2.6 X 10^ 


3.5 X 10^ 


41/50 


3.0 X 10^ 


9.8 X 10^ 


7 X 10-21 


48/50 


1.3 X 10^ 


2.7 X 10^ 


33/50 


1.6 X 10^ 


4.3 X 10^ 


5 X 10-21 


45/50 


6.7 X 10"^ 


9.5 X 10^ 


34/50 


7.7 X 10^ 


3.1 X 10^ 


3 X 10-21 


37/50 


2.5 X 10"^ 


4.1 X 10^ 


30/50 


3.1 X 10^ 


1.0 X 10^ 



Table I: Number of injections detected, mean and a for the distribution of the likelihood values for 
the GW and non-GW injections 

The major effect of this cross check is that in some cases the value of the parameter r 

0.05, which imphes 



defined in eq.(35) is above the threshold that we had chosen, rthj. 
that a direction reconstruction through the linear algebra method is simply not meaningful. 
When this happens (in most of the non-GW injections and in a few GW ones), the trigger 
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GW 
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1-d 




a 




a 




a 


10 


48/50 


0.031 


0.015 


0.037 


0.026 


0.027 


0.018 


7 


33/48 


0.048 


0.024 


0.062 


0.040 


0.040 


0.023 


5 


30/45 


0.068 


0.035 


0.089 


0.055 


0.063 


0.036 


3 


24/37 


0.11 


0.053 


0.19 


0.21 


0.15 


0.18 



Table II: Number, mean and a for the distributions of distances between injection direction and 
likelihood reconstructed (i-1), injection and determinant reconstructed (i-d), likelihood and direc- 
tion (1-d) for GW injections for different signal strengths. 





non 


-GW 




/i„,(10-2lHz-l/2) 




^J' 


a 


10 


8/41 


0.014 


0.0089 


7 


5/33 


0.050 


0.016 


5 


7/34 


0.13 


0.14 


3 


9/30 


0.13 


0.11 



Table III: Number, mean and a for the distributions of distances between likelihood and determi- 
nant reconstructed directions for non-GW injections for different signal strengths. 



is discarded. 

In Figs. T2p3 we report the distribution of distances between different direction recon- 



structions, see also tab. Ill, while fig. 14 shows graphically the result of the direction recon- 
struction with the different methods for the GW injections. As expected (see [13], [IE]), the 
precision in the determination of the arrival direction degrades at decreasing SNR. 

Finally, the application of the cut to eliminate triggers with 6 > 6thr does not change 
much the situation at this stage: we found indeed that when a direction reconstruction is 
possible, this is almost always compatible with the one found with the likelihood algorithm. 

Figure [15] shows the final detection efficiency for all the sets of injections: it should be 
noted that the efficiency for non-GW seems not to depend on the signal SNR, which points 
towards the fact that they "survived" the cuts simply because they had by chance a GW-like 
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Figure 12: Distribution of the distances between the injection direction and the direction recon- 
structed by the hkehhood method (left, filled-red), injection direction and the one reconstructed by 
the determinant method (left, empty-blue), and between likelihood and determinant reconstructed 
directions (right, empty-black), for h^ss = 10 (up) and 7 (down) •W^^^Hz^^^'^ . Also is displayed 
the distance in the direction determination for non GW-injections (right, filled-green) . 
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Figure 13: Same as Fig. 12 for hrss = 5 (up) and 3 (down) -10 ^^Hz 
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Figure 14: Each graph correspond to 50 hnearly polarized injections of a signal (25+, 25 x). The 
black circle marks the injection direction, red crosses indicate direction reconstructed with the 
likelihood method, blue dots are obtained via the determinant method. 



geometrical configuration. 



V. CONCLUSIONS 



We have set up an algorithm which exploits the multichannel capability of the sphere in 
such a way to filter spurious disturbances. We have tested our method on various sets of 
injections of GW and non-GW signals over a Gaussian background noise which reproduces 
a possible outcome of the six transducers of a spherical detector like miniGRAIL. 

We found that by crossing the two direction reconstruction methods, the false alarm 
rate due to Gaussian noise can be reduced to (9(10^^ -i- 10^'^Hz^^) and that the detection 
efficiency for non-Gaussian noise (the non-GW injections) is 5 — 10%, independently of the 
SNR. These results can be achieved while keeping the detection efficiency for GW signals in 
the range 50 — 95%. 
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Signal strength (Hz^'^) 



Figure 15: Detection efficiency of GW-injections (filled squares) and of non-GW signals (empty 
squares) as a function of the injected signal strength. The triggers surviving all the cuts are (in 
order of decreasing SNR) 48,33,30 and 21 for the GW injections, and 7,4,6 and 8 for non-GW. 



Moreover as summarized in tab. Ill, the arrival direction of the gravitational wave can 
be determined with good accuracy on the trigger satisfying our cross-check, with a precision 
which ranges from roughly 0.2 rad (~ 10°) at SNR~ 16 to rapidly improve with the signal 
strength up to 0.03 (^ 2°) at SNR~ 53. 
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Appendix A: detector parameters 



sphere mass 


JVj! 


1 1 /I Q ItTtr 


sphere radius 


-ft 


u.uozo m 


sphere temperature 


T 


0.02 K 


fundamental quadrupole modes frequencies 


^N=0,...,b 


2879,2872,2850,2741,2738 Hz 


scalar mode frequency (ioo) 




oyUz xlz 


first excited quadrupole modes frequencies 

/V2 v2 v2 V'i V'2 \ 


^N=6,...,10 


t^l^OI c; KA70 KOR'i TJrz 
00zi,00i0,t)4( Z,0ZDO,0Zi3 ( XlZ 


vector modes frequencies 


(^N=ll,...,13 


Qonn QQon qq^ji tTrj 

oyuu,ooyu,ooDi xlz 


1=3 modes frequencies 

(VZ vz vz VZ v3 vz \ 
'-'^lc'-'^ls'-'^2c'-'^2s'-'^3c'-'^3sy 


(^N=14:,...,20 


4291,4280,4247,4085,4081,4085,4081 Hz 


1=4 modes frequencies 

(Vi v4 v4 V'l V'l V'l V4 V'l V4 \ 
l-'^O )-'^lc'-'^ls'-'2c'-'^2s5^3c)-'3s'-'4c'-'4sJ 


^N=21,...,29 


c;c;nc; t^/ioi ^^aaq f^oQt; 
c^o/ii KO'iK KOA^ KO'iK vfr, 

OZ4i,0/oO,0/4i,OZoO ilZ 


sphere modes quality factors 


Qn=0...29 


2 X 10^ 


fundamental quadrupole modes 
radial displacement 


(XN=0..A 


SSQ1 1 

-z.ooyii 


scalar mode radial displacement 


"Af=5 


-o.4uyo 


first excited quadrupole modes 

icLQlcti QlbpictCcIIlcIlT; 


CKAr=6...10 


U.U I ^^ou^ 


■\7'PP'1t^t* Tn r^H pci tq H i q 1 Hi citiI q ppttipti'I' 

Vd^UVJJ- ±l±VJ^J.CO LoAJ.lClil ^J-lO^lCH^CllldlU 


'-••Af=11...13 


818002 


1=3 modes radial displacement 


«Af=14...20 


- 3.34268 


1=4 modes radial displacement 


Q;A'"=21...29 


- 3.63078 




X 


-0.3278 




X5 


-3.8043 



Table IV: Sphere features 
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transducers' azimuth 


^fc=0,l,2 


37.3773° 


transducers' azimuth 


^fe=3,4,5 


79.1876° 


transducers' longitude 


'Pk=0,l,2 


(1 + 2k) * 60° 


transducers' longitude 


¥'ifc=3,4,5 


2k * 60° 


transducers' frequencies 




2863 Hz 


transducers' frequencies 


^k=2,3 


2850 Hz 


transducers' frequencies 


^k=4,5 


2878 Hz 


transducers' masses 


'mk=o,i 


0.205 Kg 


transducers' masses 




0.153 Kg 


transducers' masses 




0.150 Kg 


transducers' quality factors 


Qk=0...6 


1 X 10^ 


gap electric fields 


Ek=0...5 


2 X lOV/m 


effective resistances 


'r'k=o,i 


0.0883 Ohm 


effective resistances 


'fk=2,3 


0.1140 Ohm 


effective resistances 


rk=4,5 


0.0872 Ohm 


total capacities 


Ck=Q,l 


1.1638x10"^ F 


total capacities 


Ck=2,3 


0.6978x10-9 F 


total capacities 


Ck=0,l 


1.1935x10-9 F 


primary inductances 


J-il n 

^^ — U...U 


0.3595 H 


secondary inductances 


T ^ 

fZ — U. . .0 


2.1x10-6 H 


auto-inductances 


K-0 5 

K — U...O 


1x10-^ H 


squid's input inductances 


fv — U. . .O 


1.7x10-6 H 


squid's self inductances 


tSQ 
^k=0...5 


8x10-11 H 


squid's mutual inductances 


^S)...5 


1x10-10 H 


squid's shunt resistances 


^f=0...5 


20 n 



Table V: Readout specifications 



